// Programa de scilab escrito 23/11/2010
// Como trabalho 4 da matéria FFI 0321 - Modelagem Matemático-Computacional
// no IFSC-USP
// _o_o_ oOo _o_o_

// A) Filtros
// Implemente, utilizando a transformada rápida de Fourrier(FFT, a função fft() do scilab) filtros passa-baixas e filtros passa-altas para sinais unidimensionais. O corte das frequências deverá ser feito de maneira abrupta, utilizando máscaras do tipo caixa, ou de maneira gradual utilizando máscaras gaussianas.

TAM=10000;
sinal = sin(100*2*%pi*(1:TAM)/8000);
SINAL = fft(sinal);
// já que a convenção é de -pi a pi,
corte = 0.5;
corte_na_fft = int((TAM/2)*corte);
resto_corte=2*(int(TAM/2) - corte_na_fft);
filtro = [ones(1,corte_na_fft), zeros(1,resto_corte), ones(1,corte_na_fft)];

SINAL_FILTRADO = filtro.*SINAL;
sinal_filtrado = real(fft(SINAL_FILTRADO,1));

function [sinal_filtrado] = filtro(sinal,freq_corte,tipo,tipo_corte)
  // filtro(sinal,0.5,"lp","gauss")
  // sinal unidimensional
  // freq_corte entre 0 e 1
  // tipo lp ou hp
  // tipo_corte "abrupto" ou "gaussiano"
  // NAO USAR SINAIS COM MENOS DE  100 AMOSTRAS
  // Para modificar a suavidade da transição gaussiana, modifique a variável barr
  TAM = length(sinal)
  corte_na_fft = int((TAM/2)*corte);
  resto_corte=(int(TAM/2) - corte_na_fft);
  
  if tipo == "lp" then // caso Low Pass
    if tipo_corte == "abrupto" then
      mascara_filtro = [ones(1,corte_na_fft), zeros(1,2*resto_corte), ones(1,corte_na_fft)];
    else
      barr = (0:.5:-.5+resto_corte/2)
      foo = exp(-barr);
      mascara_filtro = [ones(1,corte_na_fft), foo, foo($:-1:1), ones(1,corte_na_fft)];
      length(mascara_filtro)
      plot(mascara_filtro)
    end
  else // caso High Pass
    if tipo_corte == "abrupto" then
      mascara_filtro = [zeros(1,corte_na_fft), ones(1,2*resto_corte), zeros(1,corte_na_fft)];
    else
      barr=(0:.5:-.5+20)
      foo = exp(-barr)
      mascara_filtro = [zeros(1,corte_na_fft), foo($:-1:1), ones(1,2*(resto_corte)-80),foo,ones(1,corte_na_fft)];
    end
  end
  
  SINAL_FILTRADO = mascara_filtro.*SINAL;
  sinal_filtrado = real(fft(SINAL_FILTRADO,1));
endfunction

//b) Derivadas utilizando a transformada de Fourier.
//Construa uma função que, utilizando a FFT retorne a derivada de um sinal unidimensional. Como entrada é necessário o sinal a se derivar, assim como o grau da derivada, sendo este grau um número real.

SINAL_DERIVADA = %i*2*%pi*( int(-(TAM)/2):int((TAM-1)/2) ).*SINAL;
sinal_derivada=real(fft(SINAL_DERIVADA,1));
sinal_derivada2=sinal_derivada/max(sinal_derivada); // porque a fft-1 explode? TTM



//c) Aplicação sobre funções conhecidas.
//A partir de uma função conhecida (seno, cosseno) aplique os filtros implementados(passa-altas, passa-baixas e derivada) e observe quais são os resultados.
// Já está feito em a) e b)



//d) Sinais de Áudio
//A partir de um sinal de áudio, aplique os filtros passa-baixas e passa-altas e observe seus resultados.
//Dica: Utilize arquivos em formado .wav e a função do scilab waread() para ler estes dados.

[y,Fs,bits]=wavread("/usr/share/scilab/modules/sound/demos/chimes.wav")
y_lp=filtro(y,0.3,"lp","gauss")
y_hp=filtro(y,0.3,"hp","abrupto")





